Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.
In this file we have compared the 3 models HMSC,GJAM, JSDM on 21 simulated Dataset.\ We have compared the ability to recover the true interactions used for modelling this datasets, by comparing estimated correlation matrix and matrix representing true interactions.
There is a separate file, with the functions needed to fit each model and study the convergence and more parameters are presented.
## Environment filtering 5 species
data<-sim_data$EnvEvenSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
Y_cor<-cor(data$Y)
to_prec<-function(m){
n<-dim(m)[1]
Tau_n<-matrix(nrow=n, ncol=n)
for (j in 1:n) {
for (k in 1:n){
Tau_n[j, k] <- -m[j, k]/sqrt((m[j,j]*m[k,k]))
}
}
return(Tau_n)
}
#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
cols = colorRampPalette(c("dark blue","white","red"))
col2 <- colorRampPalette(c("#4393C3", "#2166AC", "#053061",
"#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
"#67001F", "#B2182B", "#D6604D", "#F4A582"))
corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")
j_metric_e5<-metrics_jsdm(me5)
#cat("Success rate ")
#j_metric_e5$success_env
cat(sprintf("Success rate: %s\n", metrics_jsdm(me5)$success_env))
## Success rate: 0.9
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =as.matrix(me5$mean$Rho*(!me5$overlap0$Rho)))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =as.matrix(me5$mean$Rho*(!me5$overlap0$Rho)))
############################################################Functions
makeSymm <- function(m) {
m[upper.tri(m)] <- t(m)[upper.tri(m)]
return(m)
}
expandSigma_rmd <- function(sigma, S){
ss <- diag(S)
ss[lower.tri(ss,diag=T)] <- sigma
ss[upper.tri(ss)] <- t(ss)[upper.tri(ss)]
ss
}
convert_to_m<-function(ar){
d <-floor((sqrt(length(ar)*8+1)-1)/2)
C <- matrix(0,d,d)
i.lwr <- which(lower.tri(C, diag = TRUE), arr.ind=TRUE)
C[i.lwr] <- ar
C<-makeSymm(C)
return(t(C))
}
metrics_gjam<-function(Rho, comp=NULL,fac=NULL,only_env=T){
if(!is.null(comp) & !is.null(fac)){
fac_comp <- fac-comp
TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac_comp[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac_comp[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac_comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac_comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
TP_comp=TP_comp+1
if(fac_comp[i,j]<0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
wrong=wrong+1
}
}
}else{
if(!is.null(comp)){
TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
TP_fac <- FN_fac<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(comp[i,j]>0 & Rho[i,j] <0)
TP_comp=TP_comp+1
if(comp[i,j]>0 & Rho[i,j] >0)
wrong=wrong+1
if(comp[i,j]>0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if(comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
if(!is.null(fac)){
TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
TP_comp <- FN_comp<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac[i,j]>0 & Rho[i,j] <0)
wrong = wrong+1
if(fac[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
}
if(only_env){
FP <- TN <- 0
TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):nrow(Rho)){
if(Rho[i,j] != 0)
FP=FP+1
if(Rho[i,j] == 0)
TN=TN+1
}
}
}
success_env<-TN/(TN+FP)
if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
"wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}
fit_gjam<-function(data, it=2500,burn=500 , name="./gjam_models/temp.rda",interact=diag(ncol(data$Y))){
#setup parameters
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
xdata<-as.data.frame(data$X[,-1])
colnames(xdata)<- c("env","env2")
ydata<-as.data.frame(data$Y)
###formula
formula<-as.formula( ~env+ env2)
ml <- list(ng = it, burnin = burn, typeNames = 'PA')
####fit
mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
save(mod_gjam1, file = name)
#summary(mod_gjam1)
hist(effectiveSize(mod_gjam1$chains$sgibbs), main="ess(beta)",lwd=2,col=gray(.6))
hist(effectiveSize(mod_gjam1$chains$bgibbs), main="ess(beta)",lwd=2,col=gray(.6))
# postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
# postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
pH<-convert_to_m(postH)
pL<-convert_to_m(postL)
R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
sgibbs<-mod_gjam1$chains$sgibbs
tau<-array(NA,dim=c(data$J,data$J,1000))
for(j in 1:1000){
ss <- expandSigma_rmd(sgibbs[j,], S = data$J)
si <- solve(ss)
tau[,,j] <- -cov2cor(si)
}
# sgibbs<-mod_gjam1$chains$sgibbs
# tau<-array(NA,dim=c(data$J,data$J,it-burn))
# for(j in (burn:it)){
# ss <- expandSigma_rmd(sgibbs[j,], S = 5)
# si <- solve(ss)
# tau[,,j-burn] <- -cov2cor(si)
# }
tau_mean<-apply(tau,c(1,2), mean)
tau_HI<-apply(tau,c(1,2),quantile,0.95)
tau_LO<-apply(tau,c(1,2),quantile,0.05)
Tau_sign<-tau_mean*(!(tau_HI>0 & tau_LO<0))
# par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
# corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Correlation cor(Y)")
# corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R")
# corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("E matrix")
# corrplot(Tau_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Tau signif")
# corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R signif")
# corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("True interactions")
R_sign
}
###this function only loads the model and return the R and T significant
load_gjam<-function(data,it=2500,burn=500,name="./gjam_models/temp.rda",interact=diag(ncol(data$Y)),autoC=F){
#setup parameters
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
gj_mod<-load_object(name)
#summary(gj_mod)
S<-ncol(data$Y)
ind <-which(array(diag(S)[lower.tri(diag(S),diag=T)])==1)
par(mfrow=c(1,3),oma = c(1, 1, 1, 1))
hist(effectiveSize(gj_mod$chains$sgibbs[,-ind]), main="sigma non diagonal elements",lwd=2,col=gray(.6))
hist(effectiveSize(gj_mod$chains$sgibbs[,ind]), main="sigma diagonal elements",lwd=2,col=gray(.6))
hist(effectiveSize(gj_mod$chains$bgibbs), main="ess(beta)",lwd=2,col=gray(.6))
if(autoC){
plot(acfplot(mcmc(gj_mod$chains$bgibbs)))
plot(acfplot(mcmc(gj_mod$chains$sgibbs)))
}
postH<-apply(gj_mod$chains$sgibbs[burn:dim(gj_mod$chains$sgibbs)[1],], 2, quantile,0.95)
postL<-apply(gj_mod$chains$sgibbs[burn:dim(gj_mod$chains$sgibbs)[1],], 2, quantile,0.05)
pH<-convert_to_m(postH)
pL<-convert_to_m(postL)
R_sign<-cov2cor(gj_mod$parameters$sigMu)*(!(pH>0 & pL<0))
sgibbs<-gj_mod$chains$sgibbs
tau<-array(NA,dim=c(data$J,data$J,it-burn))
for(j in (burn:it)){
ss <- expandSigma_rmd(sgibbs[j,], S = data$J)
si <- solve(ss)
tau[,,j-burn] <- -cov2cor(si)
}
tau_mean<-apply(tau,c(1,2), mean)
tau_HI<-apply(tau,c(1,2),quantile,0.95)
tau_LO<-apply(tau,c(1,2),quantile,0.05)
Tau_sign<-tau_mean*(!(tau_HI>0 & tau_LO<0))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(data$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(gj_mod$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(gj_mod$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("E matrix")
corrplot(Tau_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")
corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R signif")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
return(list(Rho_sign=R_sign,Tau_sign=Tau_sign))
}
######################################################################
##to setup chains parameters
data<-sim_data$EnvEvenSp5
#Rho_sign<-fit_gjam(data,2500,500,name="./gjam_models/gjam5env.rda")
mod_gjam<-load_gjam(data,name="./gjam_models/gjam5env.rda", autoC=T)
g_metric_e5<-metrics_gjam(mod_gjam$Rho_sign)
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_gjam(mod_gjam$Rho_sign)$success_env))
## Success rate: 0.5
mod_list_Rho<-list.append(mod_list_Rho,gjam=mod_gjam$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=mod_gjam$Tau)
#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="Fit",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
if (label=="Fit"){
Y_data = subset(data, select = -env)
ns<- ncol(Y_data)
np <- nrow(Y_data)
X<-scale(poly(data$env[1:np], 2))
colnames(X)<-c("env","env2")
studyDesign = data.frame(sample = as.factor(1:np))
rL = HmscRandomLevel(units = studyDesign$sample)
m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
studyDesign = studyDesign, ranLevels = list(sample = rL))
m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
save(m, file = name)
return(m)
}
if (label=="Load"){
return(load_object(name))
}
}
data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")
Convergence:
hm_conv<-function(mod, autoC=F){
codaList = convertToCodaObject(mod)
#convergence histograms
hist(effectiveSize(codaList$Beta), main="ess(beta)",lwd=2,col=gray(.6))
hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(beta)")
if(autoC){
plot(acfplot(codaList$Beta))
plot(acfplot(codaList$Omega[[1]]))
}
hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)",lwd=2,col=gray(.6))
hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(omega)")
}
hm_conv(mod=hm_mod,autoC=TRUE )
Study of interactions
metrics_hmsc<-function(Rho, comp=NULL,fac=NULL,only_env=T){
#{
if(!is.null(comp) & !is.null(fac)){
fac_comp <- fac-comp
TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac_comp[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac_comp[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac_comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac_comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
TP_comp=TP_comp+1
if(fac_comp[i,j]<0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
wrong=wrong+1
}
}
}else{
if(!is.null(comp)){
TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
TP_fac <- FN_fac<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(comp[i,j]>0 & Rho[i,j] <0)
TP_comp=TP_comp+1
if(comp[i,j]>0 & Rho[i,j] >0)
wrong=wrong+1
if(comp[i,j]>0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if(comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
# }
# }
}
if(!is.null(fac)){
TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
TP_comp <- FN_comp<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac[i,j]>0 & Rho[i,j] <0)
wrong = wrong+1
if(fac[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
}
if(only_env){
FP <- TN <- 0
TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):nrow(Rho)){
if(Rho[i,j] != 0)
FP=FP+1
if(Rho[i,j] == 0)
TN=TN+1
}
}
}
success_env<-TN/(TN+FP)
if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
"wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}
hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
getOmega = function(a,r=1)
return(crossprod(a$Lambda[[r]]))
ns<-mod$ns
postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
postOmega<-abind(postOmega1,postOmega2,along=3)
postOmegaMean = apply(postOmega,c(1,2),mean)
postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
postR<-array(dim=c(ns,ns,nchains*nsamples))
postT<-array(dim=c(ns,ns,nchains*nsamples))
for(i in 1:dim(postOmega)[3]){
postR[,,i]<-stats::cov2cor(postOmega[,,i])
postT[,,i]<- -stats::cov2cor(solve(postOmega[,,i]+diag(ns)))
}
postRMean = apply(postR,c(1,2),mean)
postRUp=apply(postR,c(1,2),quantile,0.95)
postRLo=apply(postR,c(1,2),quantile,0.05)
postTMean = apply(postT,c(1,2),mean)
postTUp=apply(postT,c(1,2),quantile,0.95)
postTLo=apply(postT,c(1,2),quantile,0.05)
Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
Toplot_T<-postTMean*(!(postTUp>0 & postTLo<0))
# Omegacor<- computeAssociations(m)
# supportLevel<- 0.95
# toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
# corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(hm_mod$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(postRMean, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(Toplot_R, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Plot only non zero value")
corrplot(Toplot_T, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation matrix")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
return(list(Rho_sign=Toplot_R,Tau_sign=Toplot_T))
}
hm_mod_R<-hm_inter(hm_mod,nchains=2,nsamples = 3000)
h_metric_e5<-metrics_hmsc(hm_mod_R$Rho_sign)
# cat("Success rate ")
# h_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_hmsc(hm_mod_R$Rho_sign)$success_env))
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R JSDM ")
corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R GJAM")
corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R HMSC")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(5))
ALL4<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
#corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
#title("Partial correlation JSDM ")
corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation GJAM")
corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation HMSC")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(5))
me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)
## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
jsdm_conv(me10)
## Maximum Rhat value for Beta: 1.5602
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.7953
me10$mcmc.info[1:7]
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
j_metric_e10<-metrics_jsdm(me10)
cat("Success rate ")
## Success rate
j_metric_e10$success_env
## [1] 0.9333333
data<-sim_data$EnvEvenSp10
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
#plot_cor_jsdm(me10,data$Y)
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =me10$mean$Rho*(!me10$overlap0$Rho))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =me10$mean$Rho*(!me10$overlap0$Rho))
data<-sim_data$EnvEvenSp10
#Rho_sign<-fit_gjam(data,3000,500,"./gjam_models/gjam10env.rda")
gjam_mod<-load_gjam(data,3000,500,name="./gjam_models/gjam10env.rda")
#gje10<-load_object("./gjam_models/gjam10env.rda")
g_metric_e10<-metrics_gjam(gjam_mod$Rho_sign)
cat("Success rate ")
## Success rate
g_metric_e10$success_env
## [1] 0.3111111
#to check posterior density of s in Sigma
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)
data<-sim_data$EnvEvenSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10env.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod)
h_metric_e10<-metrics_hmsc(hm_mod_R$Rho_sign)
# cat("Success rate ")
# h_metric_e10$success_env
# h_metric_e10$success_fac
cat(" ")
cat(sprintf("Success rate: %s\n", h_metric_e10$success_env))
## Success rate: 1
cat(sprintf("Success rate: %s\n", h_metric_e10$success_fac))
mod_list_Rho<-list.append(mod_list_Rho, hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(10))
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(10))
me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)
## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
me20$mcmc.info[1:7]
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
jsdm_conv(me20)
## Maximum Rhat value for Beta: 1.2174
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3236
j_metric_e20<-metrics_jsdm(me20)
cat("Success rate ")
## Success rate
j_metric_e20$success_env
## [1] 0.9894737
data<-sim_data$EnvEvenSp20
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
plot_cor_jsdm(me20,data$Y)
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =me20$mean$Rho*(!me20$overlap0$Rho))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =me20$mean$Rho*(!me20$overlap0$Rho))
data<-sim_data$EnvEvenSp20
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
gjam_mod<-load_gjam(data,5000,500,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")
g_metric_e20<-metrics_gjam(gjam_mod$Rho_sign)
cat("Success rate ")
## Success rate
#g_metric_e20$success_env
cat(" ")
cat(sprintf("Success rate env: %s\n", round(metrics_gjam(gjam_mod$Rho_sign)$success_env,3)))
## Success rate env: 0.358
#to check posterior density of s in Sigma
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
#
# data<-sim_data$EnvEvenSp20
#
# #fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
#
# #load_gjam(data,name="./gjam_models/gjam20env.rda")
# #gje20<-load_object("./gjam_models/gjam20env.rda")
#
#
# #to check posterior density of s in Sigma
# #gje20<-load_object("./gjam_models/gjam20env.rda")
# #plot(density(gje20$chains$sgibbs[,4]))
# data <- list(
# Y = subset(data, select = -env),
# X = cbind(1, scale(poly(data$env, 2))),
# covx = cov(cbind(1, scale(poly(data$env, 2)))),
# K = 3,
# J = ncol(data) - 1,
# n = nrow(data),
# I = diag(ncol(data) - 1),
# df = ncol(data)
# )
# xdata<-as.data.frame(data$X[,-1])
# colnames(xdata)<- c("env","env2")
# ydata<-as.data.frame(data$Y)
# ###formula
# rl <- list(r = 8, N = 20)
# formula<-as.formula( ~env+ env2)
# ml <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
# ####fit
#
#
# mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
# save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
# summary(mod_gjam1)
#
# Tau <- solve(mod_gjam1$parameters$sigMu)
# Tau_n = to_prec(Tau)
#
# #postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
# #postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
#
#
# #pH<-convert_to_m(postH)
# #pL<-convert_to_m(postL)
#
# #R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
#
# g_metric_e20_dim_red<-metrics_gjam(Rho_sign)
# cat("Success rate ")
# g_metric_e20_dim_red$success_env
#
# par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
# corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Correlation cor(Y)")
# corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R")
# corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("E matrix")
# # corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# # title("Tau")
# # corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# # title("R signif")
# corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("True interactions")
#
data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod)
h_metric_e20<-metrics_hmsc(hm_mod_R$Rho_sign)
cat("Success score ")
## Success score
h_metric_e20$success_env
## [1] 1
cat(" ")
cat(sprintf("Success rate non-interacting: %s\n", round(metrics_hmsc(hm_mod_R$Rho_sign)$success_env,3)))
## Success rate non-interacting: 1
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(20))
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(20))
mfd5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mfd5)
## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
#mfd5$Rhat
jsdm_conv(mfd5)
## Maximum Rhat value for Beta: 1.3081
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1421
mfd5$mcmc.info[1:7]
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
j_metric_facDense5<-metrics_jsdm(mfd5,fac = fac_inter[[4]],only_env = F)
cat("Success rate ")
## Success rate
j_metric_facDense5$success_env
## [1] 1
j_metric_facDense5$success_fac
## [1] 0
data<-sim_data$FacDenseSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
#plot_cor_jsdm(mfd5,data$Y,fac_inter[[4]])
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =mfd5$mean$Rho*(!mfd5$overlap0$Rho))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =mfd5$mean$Rho*(!mfd5$overlap0$Rho))
data<-sim_data$FacDenseSp5
#Rho_sign<-fit_gjam(data,10000,2000,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
gjam_mod<-load_gjam(data,10000,2000,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
g_metric_facDense5<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[4]], only_env = F)
#cat("Success rate ")
#g_metric_facDense5$success_env
#g_metric_facDense5$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0.5
mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)
data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod, interact = fac_inter[[4]])
## $Rho_sign
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1
##
## $Tau_sign
## [,1] [,2] [,3] [,4] [,5]
## [1,] -1 0 0 0 0
## [2,] 0 -1 0 0 0
## [3,] 0 0 -1 0 0
## [4,] 0 0 0 -1 0
## [5,] 0 0 0 0 -1
hm_mod_R<-hm_inter(hm_mod)
h_metric_facDense5<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[4]],only_env = F)
#cat("Success rate ")
# h_metric_facDense5$success_env
# h_metric_facDense5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[4]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[4]])
## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 0.925
## [1] 0
data<-sim_data$FacDenseSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
gjam_mod<- load_gjam(data,5000,500,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
g_metric_facDense10<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[5]], only_env = F)
# cat("Success rate ")
# g_metric_facDense10$success_env
# g_metric_facDense10$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense10$success_env,3)))
## Success rate for non-interacting: 0.625
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense10$success_fac,3)))
## Success rate for facilitation: 0.8
mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)
## Success rate for non-interacting: 1
## Success rate for facilitation: 0
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[5]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[5]])
## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9833333
## [1] 0
data<-sim_data$FacDenseSp20
#Rho_sign<-fit_gjam(data,10000,1500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
gjam_mod<-load_gjam(data,10000,1500,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
g_metric_facDense20<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[6]], only_env = F)
# cat("Success rate ")
# g_metric_facDense20$success_env
# g_metric_facDense20$success_fac
#
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 0.406
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0.5
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[6]])
h_metric_facDense20<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[6]],only_env = F)
# cat("Success rate ")
# h_metric_facDense20$success_env
# h_metric_facDense20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[6]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[6]])
## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.889
## Success rate for facilitation: 1
data<-sim_data$FacSparseSp5
#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam5fs.rda",interact=fac_inter[[7]])
gjam_mod<-load_gjam(data,2500,500,name="./gjam_models/gjam5fs.rda", interact=fac_inter[[7]])
#gjfd5<-load_object("./gjam_models/gjam5fs.rda")
g_metric_facSparse5<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[7]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse5$success_env
# g_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacSparseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fs.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[7]])
h_metric_facSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[7]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[7]])
## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.93
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp10
#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam10fs.rda",interact=fac_inter[[8]])
gjam_mod<-load_gjam(data,name="./gjam_models/gjam10fs.rda", interact=fac_inter[[8]])
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
g_metric_facSparse10<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[8]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse10$success_env
# g_metric_facSparse10$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 0.326
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacSparseSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10fs.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[8]])
h_metric_facSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse10$success_env
# h_metric_facSparse10$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[8]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[8]])
## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.978
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp20
#fit_gjam(data,10000,1500,"./gjam_models/gjam20fs.rda",interact=fac_inter[[9]])
gjam_mod <-load_gjam(data,10000,1500,name="./gjam_models/gjam20fs.rda", interact=fac_inter[[9]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
g_metric_facSparse20<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[9]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse20$success_env
# g_metric_facSparse20$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 0.446
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0.25
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacSparseSp20
hm_mod<-fit_hmsc(data,nsamples=3000, nchains=2,"Load",name="./HMmodels/hm20fs.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = fac_inter[[9]])
h_metric_facSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[9]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse20$success_env
# h_metric_facSparse20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[9]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[9]])
## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## Successful convergence based on Rhat values (all < 1.1).
## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 1
## [1] 1
## Success rate for non-interacting: 1
## Success rate for facilitation: 0
data<-sim_data$CompDenseSp5
#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam5cmpd.rda",interact=comp_inter[[10]])
gjam_mod<-load_gjam(data,10000,1000,name="./gjam_models/gjam5cmpd.rda", interact=(-1)*comp_inter[[10]])
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
g_metric_compDense5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[10]], only_env = F)
# cat("Success rate ")
# g_metric_compDense5$success_env
# g_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 0.625
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$CompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hm5cmpd.rda" )
#hm_mod<-fit_hmsc(data,"Fit",nsamples=3000, nchains=2 )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[10]])
h_metric_compDense5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[10]],only_env = F)
# cat("Success rate ")
# h_metric_compDense5$success_env
# h_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[10]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[10]])
## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.875
## Success rate for facilitation: 1
data<-sim_data$CompDenseSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10cmpd.rda",interact=comp_inter[[11]])
gjam_mod<-load_gjam(data,5000,500,name="./gjam_models/gjam10cmpd.rda", interact= (-1)*comp_inter[[11]])
#gjfd5<-load_object("./gjam_models/gjam10cmpd.rda")
g_metric_compDense10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[11]], only_env = F)
# cat("Success rate ")
# g_metric_compDense10$success_env
# g_metric_compDense10$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.325
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$CompDenseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm10cmpd.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod,nsamples=5000, nchains=2, interact = (-1)*comp_inter[[11]])
h_metric_compDense10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[11]],only_env = F)
# cat("Success rate ")
# h_metric_compDense10$success_env
# h_metric_compDense10$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.875
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[11]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[11]])
## Summary for model '/tmp/RtmpKix4lZ/file40e96843124a'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2061.424 minutes at time 2019-04-16 03:39:17.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.9065
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.6944
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.956
## Success rate for competition: 0.7
data<-sim_data$CompDenseSp20
#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam20cmpd.rda",interact= (-1)*comp_inter[[12]])
gjam_mod<-load_gjam(data,10000,1000,name="./gjam_models/gjam20cmpd.rda", interact= (-1)*comp_inter[[12]])
#gjfd5<-load_object("./gjam_models/gjam20cmpd.rda")
g_metric_compDense20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[12]], only_env = F)
# cat("Success rate ")
# g_metric_compDense20$success_env
# g_metric_compDense20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.394
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.9
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$CompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm20cmpd.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=5000, nchains=2,interact = (-1)*comp_inter[[12]])
h_metric_compDense20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[12]],only_env = F)
# cat("Success rate ")
# h_metric_compDense20$success_env
# h_metric_compDense20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.967
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.4
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[12]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[12]])
## Summary for model '/tmp/RtmpKix4lZ/file40e97e225117'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.251 minutes at time 2019-04-17 14:00:46.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1876
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2038
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 1
## Success rate for competition: 1
data<-sim_data$CompSparseSp5
#Rho_sign<-fit_gjam(data,5000,2000,"./gjam_models/gjam5cmps.rda",interact= (-1)*comp_inter[[13]])
gjam_mod<-load_gjam(data,5000,2000,name="./gjam_models/gjam5cmps.rda", interact= (-1)*comp_inter[[13]])
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
g_metric_compSparse5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[13]], only_env = F)
cat("Success rate ")
## Success rate
g_metric_compSparse5$success_env
## [1] 0.5555556
g_metric_compSparse5$success_comp
## [1] 1
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse5$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$CompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm5cmps.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[13]])
h_metric_compSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[13]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse5$success_comp,3)))
## Success rate for competition: 1
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[13]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[13]])
## Summary for model '/tmp/RtmpKix4lZ/file40e978d641ed'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 685.719 minutes at time 2019-04-17 14:30:01.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.116
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1753
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9534884
## [1] 1
## Success rate for non-interacting: 0.953
## Success rate for competition: 1
data<-sim_data$CompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam10cmps.rda",interact= (-1)*comp_inter[[14]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjam10cmps.rda", interact= (-1)*comp_inter[[14]])
g_metric_compSparse10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[14]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse10$success_env
# g_metric_compSparse10$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse10$success_env,3)))
## Success rate for non-interacting: 0.326
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$CompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm10cmps.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[14]])
h_metric_compSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[14]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse10$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse10$success_comp,3)))
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[14]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[14]])
## Summary for model '/tmp/RtmpKix4lZ/file40e9379ee963'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2084.707 minutes at time 2019-04-18 01:55:46.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2652
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3883
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.968
## Success rate for competition: 0
data<-sim_data$CompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam20cmps.rda",interact= (-1)*comp_inter[[15]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjam20cmps.rda", interact= (-1)*comp_inter[[15]])
g_metric_compSparse20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[15]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse20$success_env
# g_metric_compSparse20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 0.269
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0.75
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$CompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm20cmps.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[15]])
h_metric_compSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[15]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[15]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[15]])
## Summary for model '/tmp/RtmpKix4lZ/file40e93b49ad29'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 28.267 minutes at time 2019-04-19 12:40:31.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1667
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2086
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.667
## Success rate for competition: 1
## Success rate for facilitation: 0.5
data<-sim_data$FacCompDenseSp5
#Rho_sign<-fit_gjam(data,30000,15000,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
gjam_mod<-load_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign3<-fit_gjam(data,5000,1000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign4<-fit_gjam(data,5000,2500,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign5<-fit_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmp_facd.rda", (-1)*comp_inter[[16]]+fac_inter[[16]])
g_metric_FacCompDense5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[16]], fac=fac_inter[[16]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense5$success_fac,3)))
## Success rate for facilitation: 1
# par(mfrow=c(2,3))
#corrplot((Rho_sign+Rho_sign2)/2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign3, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign4, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign5, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacCompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hmcmp_facd5.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[16]]+ fac_inter[[16]])
h_metric_FacCompDense5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[16]],fac= fac_inter[[16]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense5$success_comp,3)))
## Success rate for competition: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense5$success_fac,3)))
## Success rate for facilitation: 0.5
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])
## Summary for model '/tmp/RtmpKix4lZ/file40e93b173410'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 688.628 minutes at time 2019-04-19 13:08:47.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.3458
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3195
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.914
## Success rate for competition: 0
## Success rate for facilitation: 0
data<-sim_data$FacCompDenseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd10.rda",interact= (-1)*comp_inter[[17]]+fac_inter[[17]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facd10.rda", interact= (-1)*comp_inter[[17]]+fac_inter[[17]])
g_metric_FacCompDense10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[17]], fac=fac_inter[[17]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense10$success_env,3)))
## Success rate for non-interacting: 0.286
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense10$success_comp,3)))
## Success rate for competition: 0.6
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense10$success_fac,3)))
## Success rate for facilitation: 1
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd10.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[17]] + fac_inter[[17]])
h_metric_FacCompDense10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[17]],fac=fac_inter[[17]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense10$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])
## Summary for model '/tmp/RtmpKix4lZ/file40e95eadc1e3'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2073.109 minutes at time 2019-04-20 00:37:26.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.3008
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2755
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.965
## Success rate for competition: 0.333
## Success rate for facilitation: 0.333
data<-sim_data$FacCompDenseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd20.rda",interact= (-1)*comp_inter[[18]]+fac_inter[[18]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facd20.rda", interact= (-1)*comp_inter[[18]]+fac_inter[[18]])
g_metric_FacCompDense20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[18]], fac=fac_inter[[18]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense20$success_env,3)))
## Success rate for non-interacting: 0.366
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense20$success_comp,3)))
## Success rate for competition: 0.667
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense20$success_fac,3)))
## Success rate for facilitation: 0.429
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacCompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd20.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[18]] + fac_inter[[18]])
h_metric_FacCompDense20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[18]],fac=fac_inter[[18]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense20$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense20$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])
## Summary for model '/tmp/RtmpKix4lZ/file40e9896a5fb'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.285 minutes at time 2019-04-21 11:10:35.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.3771
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3486
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.875
## Success rate for competition: 1
## Success rate for facilitation: 1
data<-sim_data$FacCompSparseSp5
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs5.rda",interact= (-1)*comp_inter[[19]]+ fac_inter[[19]] )
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs5.rda", interact= (-1)*comp_inter[[19]]+ fac_inter[[19]])
g_metric_FacCompSparse5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[15]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse5$success_comp,3)))
## Success rate for competition: 1
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacCompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs5.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[19]]+ fac_inter[[19]])
h_metric_FacCompSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[19]],fac= fac_inter[[19]] ,only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse5$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])
## Summary for model '/tmp/RtmpKix4lZ/file40e95f910eba'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 679.081 minutes at time 2019-04-21 11:39:53.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 2.1146
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 3.3048
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.951
## Success rate for competition: 0
## Success rate for facilitation: 0
data<-sim_data$FacCompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs10.rda",interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs10.rda", interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
g_metric_FacCompSparse10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse10$success_env,3)))
## Success rate for non-interacting: 0.341
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse10$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparse10$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacCompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs10.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[20]] +fac_inter[[20]])
h_metric_FacCompSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse10$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])
## Summary for model '/tmp/RtmpKix4lZ/file40e978c52358'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2083.296 minutes at time 2019-04-21 22:58:59.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2707
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3875
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.984
## Success rate for competition: 0
## Success rate for facilitation: 0
data<-sim_data$FacCompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs20.rda",interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs20.rda", interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
g_metric_FacCompSparse20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[21]], fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse20$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse20$success_comp,3)))
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparse20$success_fac,3)))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign)
data<-sim_data$FacCompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs20.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[21]] +fac_inter[[21]])
h_metric_FacCompSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[21]],fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse20$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse20$success_comp,3)))
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse20$success_fac,3)))
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])